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PRINCIPAL NOTATION 


SYMBOLS 

Unless specified otherwise, all variables are nondimensional. 
Symbol Definition 


A, B, C 
A', B', C' 

E, F 

A A 

E, F 
E T 

E K , F y 

A A 

E y, F y 

E L, 

A A 

Ek 2 , f , 2 
K 

H, H„ 

H, Hy 

ij 

J 

k 

k„ k, 

L, 

K 

N» N 2 

P 

Pr r 

Pr„ Pr, 


Speed of sound. 

Coefficient submatrices in block tridiagonal system of equations. 

Coefficient submatrices for boundary conditions. 

Specific heats at constant pressure and volume. 

Inviscid flux vectors in the Cartesian or cylindrical coordinate form of the govern- 
ing equations. 

Inviscid flux vectors in the computational coordinate form of the governing 
equations. 

Total energy per unit volume. 

Viscous flux vectors in the Cartesian or cylindrical coordinate form of the govern- 
ing equations. 

Viscous flux vectors in the computational coordinate form of the governing 
equations. 

Non-cross derivative viscous flux vectors in the computational coordinate form of 
the governing equations. 

Cross derivative viscous flux vectors in the computational coordinate form of the 
governing equations. 

Stagnation enthalpy per unit mass. 

Non-derivative inviscid and viscous terms in the Cartesian coordinate form of the 
governing equations for axisymmetric flow. 

Non-derivative inviscid and viscous terms in the computational coordinate form 
of the governing equations for axisymmetric flow. 

Grid indices in the £ and y\ directions. 

Jacobian matrix of the generalized grid transformation. 

Effective thermal conductivity coefficient. 

Laminar and turbulent thermal conductivity coefficient. 

Dimensional reference length. 

Number of governing equations being solved. 

Number of grid points in the £ and rj directions. 

Static pressure. 

Reference Prandtl number. 

Laminar and turbulent Prandtl number. 

Heat fluxes in the cylindrical x and r directions. 
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Definition 


Symbol 

Q 

Q 

R 

Re, 

S 

S' 

t 

T 

u, v 

w, V, w 

X, r 

x, y 


y 

6 

A, V 

4 2 \ 4 4) 


ef\ ef\ etc. 

01 , 02 . ^ 

K 2i *4 

A 

h K 

M 

M/> Mr 
V 

£, >i 
p 


Tjt*» T^, etc. 
* 


Heat fluxes in the Cartesian jc and j; directions. 

Vector of dependent variables in the Cartesian or cylindrical coordinate form of the 
governing equations. 

Vector of dependent variables in the computational coordinate form of the gov- 
erning equations. 

Gas constant. 

Reference Reynolds number. 

Source term subvector in block tridiagonal system of equations. 

Source term subvector for boundary conditions. 

Physical time. 

Static temperature. 

Velocities in the Cartesian x and y directions. 

Velocities in the cylindrical x x r, and swirl directions. 

Cylindrical axial and radial coordinates. 

Cartesian coordinates. 

Centering parameter in differencing formula for spatial first derivatives. 

Ratio of specific heats, c p /c v . 

Difference operator. 

First-order forward and backward difference operators. 

Second- and fourth-order explicit artificial viscosity coefficients in constant coeffi- 
cient model. 

Implicit artificial viscosity coefficient. 

Second- and fourth-order artificial viscosity coefficients in nonlinear coefficient 
model. 

Parameters determining type of time differencing used. 

Constants in nonlinear coefficient artificial viscosity model. 

Effective second coefficient of viscosity. 

Laminar and turbulent second coefficient of viscosity. 

Effective viscosity coefficient. 

Laminar and turbulent viscosity coefficient. 

Laminar kinematic viscosity. 

Computational coordinate directions. 

Static density. 

Pressure gradient scaling parameter in nonlinear coefficient artificial viscosity 
model. 

Computational time. 

Elements of shear stress tensor. 

Spectral radius in nonlinear coefficient artificial viscosity model. 


4 Principal Notation 


PROTEUS 2-D Analysis Description 



SUBSCRIPTS 


Subscript 

r 

t 

x> r 
L n 

T 

SUPERSCRIPT? 

Superscript 

n 

+ 


Definition 

Denotes grid location in ^ and r\ directions. 

Denotes dimensional reference condition. 

Denotes differentiation with respect to physical time. 

Denotes differentiation with respect to cylindrical coordinate directions. 
Denotes differentiation with respect to Cartesian coordinate directions. 
Denotes differentiation with respect to computational coordinate directions. 
Denotes differentiation with respect to computational time. 


Definition 
Denotes time level. 

Denotes solution after first ADI sweep. 
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SUMMARY 

A new computer code, called PROTEUS, has been developed to solve the two-dimensional planar or 
axisvmmetric, Reynolds-averaged, unsteady compressible Navier-Stokes equations in strong conservation 
law form The objective in this effort has been to develop a code for aerospace propulsion applications that 
is easy to use and easy to modify. Code readability, modularity, and documentation have been emphasized. 

The governing equations are written in Cartesian coordinates and transformed into generalized 
nonorthogonal body-fitted coordinates. They are solved by marching m time using a fully-coupled 
altemating-direction-implicit solution procedure with generalized first- or second-order time differencing. 
The boundary conditions are also treated implicitly, and may be steady or unsteady. Spatially periodic 
boundary conditions are also available. All terms, including the diffusion terms are linearized using 
second-order Taylor series expansions. Turbulence is modeled using an algebraic eddy viscosity model. 

The program contains many operating options. The governing equations may be solved for two- 
dimensional planar flow, or axisymmetric flow with or without swirl. The thin-layer or Euler equa ions 
may be solved as subsets of the Navier-Stokes equations. The energy equation may be eliminated by the 
assumption of constant total enthalpy. Explicit and implicit artificial viscosity may be used to damp pre- 
and post-shock oscillations in supersonic flow and to minimize odd-even decoupling caused by central 
spatial differencing of the convective terms in high Reynolds number flow. Several tune step options are 
available for convergence acceleration, including a locally variable time step and global time step cycling. 
Simple Cartesian or polar grids may be generated internally by the program. More complex geometnes 
require an externally generated computational coordinate system. 

The documentation is divided into three volumes. Volume 1, the current volume is the Analysis De- 
scription, and presents the equations and solution procedure used m PROTEUS. It describes in detail the 
governing equations, the turbulence model, the linearization of the equations and boundary conditions the 
time and space differencing formulas, the ADI solution procedure, and the artificial viscosity models. Vol- 
ume 2 is the User's Guide, and contains information needed to run the program. It desenbes the program s 
general features, the input and output, the procedure for setting up initial conditions, the computer resource 
requirements, the diagnostic messages that may be generated, the job control language used to run e 
program, and several test cases. Volume 3 is the Programmer's Reference, and contains detailed informa- 
tion useful when modifying the program. It describes the program structure, the Fortran vanables stored 
in common blocks, and the details of each subprogram. 
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1.0 INTRODUCTION 


Much of the effort in applied computational fluid dynamics consists of modifying an existing program 
for whatever geometries and flow regimes are of current interest to the researcher Unfortunately, near y 
all of the available nonproprietary programs were started as research projects with the emphasis on dem- 
onstrating the numerical algorithm rather than ease of use or ease of modification I he developers u^ally 
intend to clean up and formally document the program, but the immediate need to extend it to new ge- 
ometries and flow regimes takes precedence. 

The result is often a haphazard collection of poorly written code without any consistent structure. An 
extensively modified program may not even perform as expected under certain combinations of operating 
options. Each new user must invest considerable tune and effort in attempting to understand the underlying 
structure of the program if intending do anything more than run standard test cases with it. The user s 
subsequent modifications further obscure the program structure and therefore make it even more difficult 
for others to understand. 

The PROTEUS two-dimensional Navier-Stokes computer program is a user-oriented and easily- 
modifiable flow analysis program for aerospace propulsion applications. Readability, modularity, and 
documentation were primary objectives during its development. The entire program was specified, de- 
signed, and implemented in a controlled, systematic manner. Strict programming standards were enforced 
by immediate peer review of code modules; Kemighan and Plauger (1978) provided many useful ideas about 
consistent programming style. Every subroutine contains an extensive comment section describing the 
inpu.^31 oi.pu. variables, and calling sequence of the subroutine. Wrth ,us. two clearly- 
defined exceptions, the entire program is wntten in ANSI standard Fortran 77 to enhance portability. / 
master version of the program is maintained and periodically updated with corrections, as well as extensions 
of general interest (e.g., turbulence models.) 

The PROTEUS program solves the unsteady, compressible, Reynolds-averaged Navier-Stokes 
equations in strong conservation law form. The governing equations are written in Cartesian coordinates 
and transformed into generalized nonorthogonal body-fitted coordinates. They are solved by marching in 
time using a fully-coupled altemating-direction-implicit (ADI) scheme with generalized time and space dif- 
ferencing (Briley and McDonald, 1977; Beam and Warming, 1978). The current turbulence model is based 
upon the algebraic eddy- viscosity model of Baldwin and Urniax (1978). All tenns, including the diffusion 
terms, are linearized using second-order Taylor series expansions. The boundary conditions are treated 
implicitly, and may be steady or unsteady. Spatially periodic boundary conditions are also available. 

The program contains many operating options. The governing equations may be solved for two- 
dimensional planar flow, or axisymmctric flow with or without swirl. The thin-layer or Euler equations 
may be solved as subsets of the Navier-Stokes equations. The energy equation may' be eliminated by the 
assumption of constant total enthalpy. Explicit and implicit artificial viscosity may be used to damp pre- 
and post-shock oscillations in supersonic flow and to minimize odd-even decoupling caused by central 
spatial differencing of the convective terms in high Reynolds number flow. Several tune step options are 
available for convergence acceleration, including a locally variable tune step and global tune step cycling. 
Simple grids may be generated internally by the program; more complex geometnes requue external grid 
generation, such as that developed by Chen and Schwab (1988). 

The documentation is divided into three volumes. Volume 1, the current volume, is the Analysis De- 
scription, and presents the equations and solution procedure used m PROTEUS. It describes in detail the 
governing equations, the turbulence model, the linearization of the equations and boundary conditions the 
time and space differencing formulas, the ADI solution procedure, and the artificial viscosity models. Vol- 
ume 2 is the User's Guide, and contains information needed to run the program. It desenbes the program s 
genera] features, the input and output, the procedure for setting up initial conditions, the computer resource 
requirements, the diagnostic messages that may be generated, the job control language used to run the 
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program and several test cases. Volume 3 is the Programmer's Reference, and contains detailed informa- 

J* " “f ful when modifying the program It describes the program structure, the Fortran variables stored 
in common blocks, and the details of each subprogram. 

The authors would like to acknowledge the significant contributions made by three co-workers in the 
development of the PROTEUS program. Simon Chen did the original coding of the Baldwin-Lomax tur- 
SS* £° de [’ ^nd consulted in the implementation of the nonlinear coefficient artificial viscosity model 
William Kumk developed the original codmg for computing the metrics of the generalized nonorthogonal 

SErioSSZSy K. madC m “ y debug8ing and veHr ' ca,ion Particularly for spatiaUy 
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2.0 GOVERNING EQUATIONS 


2.1 GOVERNING EQUATIONS IN CARTESIAN COORDINATES 

The basic governing equations are the two-dimensional compressible Navier-Stokcs equations. These 
equations may be found in several standard references (e.g., Hughes and Gaylord, 1964; Schlichting, 1968; 
White, 1974; Anderson, Tannehill, and Pletcher, 1984.) In Cartesian coordinates, the two-dimensional 
planar equations 1 can be written in strong conservation law form using vector notation as 

iQ ^ = + (2.1) 

dt dx dy dx dy 


where 


Q = [p P« P v e t\ T 


pu 

pu 2 +p 

puv 

( E t + p)u 


pv 

puv 
pv 2 +p 
(E r + p)v 


(2.2a) 


(2.2b) 


(2.2c) 


E>/ = 


Re r 


(2.2d) 


l xy 


ur xx + vr xy p r 


4x\ 


Fi/ = 


Re r 


‘xy 


iyy 


UT xy VT yy p r $y 


(2-2e) 


Equation (2.1) thus represents, in order, the continuity, x-momentum, j;-momentum, and energy equations, 
with dependent variables p y pu , pv, and E T . 


1 PROTEUS can be used for both two-dimensional planar or axisymmetric flow. However, the axisymmetric 
equations have some additional terms that complicate the analysis somewhat. For the sake of clarity, the main 
body of this report describes the two-dimensional planar analysis, and the axisymmetric analysis is described in 
Appendix B. 
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The shear stresses and heat fluxes are given by 



In these equations, t represents time; x and y represent the Cartesian coordinate directions; u and v are 
the velocities in the x and y directions; p, /?, and T are the static density, pressure, and temperature; E r is 
the total energy' per unit volume; and p, A, and k are the coefficient of viscosity, second coefficient of 
viscosity, and coefficient of thermal conductivity. 

All of the above equations have been nondimensionalized using appropriate normalizing conditions. 
Lengths have been nondimensionalized by L r , velocities by u n density by p n temperature by T n viscosity 
by thermal conductivity by k n pressure and total energy by p r u } } , and time by LJu r The reference 
Reynolds and Prandtl numbers are thus defined as Re r — p r u r L r jp r and Pr f — p r u?/k r T r . 2 

Turbulence is modeled using the Boussinesq approach (Schlichting, 1968). The equations presented in 
this section are thus used for both laminar and turbulent flow. For turbulent flow they represent the 
Reynolds time-averaged form of the Navier-Stokes equations, with density fluctuations neglected. They 
may also be interpreted as the Favre or mass-weighted time-averaged form of the equations. With Favre 
time averaging, however, the velocities and thermal variables represent mass-averaged quantities defined by 
u — pujp, etc., where the overbar represents a conventional Reynolds time-averaged quantity. Details on 
Reynolds and Favre time- averaging procedures may be found in Cebeci and Smith (1974), and in Anderson, 
Tannehill, and Pletcher (1984). In either case, p, A t and k represent effective coefficients. For example, in 
turbulent flow \i = \i t 4- p n where and are the laminar and turbulent viscosity coefficients, and p t comes 
from some appropriate turbulence model. The model currently used in the PROTEUS code is the algebraic 
eddy viscosity model of Baldwin and Lomax (1978), implemented as described in Section 3.0. 

2.2 EQUATION OF STATE 

In addition to the equations presented above, an equation of state is required to relate pressure to the 
dependent variables. Any appropriate equation, or even table, could be used. The equation currently built 
into the PROTEUS code is the equation of state for thermally perfect gases, p = pRT, where R is the gas 
constant. For calorically perfect gases, this can be rewritten as 

p = (y - l)[f r - -i- p(u 2 + v 2 )] (2.4) 

where y is the ratio of specific heats, c p lc v . Here the gas constant and specific heats have been 
nondimensionalized by w r 2 /T r . 

If the flow is such that we can assume a perfect gas with constant stagnation enthalpy, the energy 
equation may be eliminated. This assumption is reasonable, for example, in inviscid regions, and in 


2 Note that this Prandtl number does not have a physically meaningful value, but is merely defined by a combination 
of the normalizing conditions for c pr and k that appear when the equations are nondimensionalized. 
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adiabatic wall boundary layers if the Prandtl number is near 1 (Briley and McDonald, 1977). The stag- 
nation enthalpy is defined as 

h T =c p T+±.(u 2 + v 2 ) (2.5) 

Here the stagnation enthalpy is nondimensionalized by uj. The temperature is thus 

T = -jr- [h T - y (u 2 + v 2 ) J (2.6) 


and the equation of state becomes 

p = y - y~ - p [v - y (“ 2 + y2 )] ( 2 - 7 ) 

This equation of state does not require the total energy E T , and the energy equation need not be solved. 
The total energy may be computed from 

E T = ph T — p (2.8) 


2.3 GENERALIZED GRID TRANSFORMATION 


Because the governing equations in the previous section are written in Cartesian coordinates, they are 
not well suited for general geometric configurations. For most applications a body-fitted coordinate system 
is desired. This greatly simplifies the application of boundary conditions and the bookkeeping in the nu- 
merical method used to solve the equations. The following generalized grid transformation, which can be 
orthogonal or nonorthogonal, is therefore used to transform the governing equations from physical {x,y, t) 
coordinates to rectangular orthogonal computational (£, tj, t) coordinates. 

t = Z(x,y. 0 

= 0 ( 2 - 9 ) 

T = t 

In PROTEUS, the spatial computational domain is square, with £ and rj each running from 0 to 1. Using 
the chain rule for partial differentiation, the derivatives in the Cartesian form of the governing equations can 
be replaced using the following expressions. 


_d_ , _d_ 

dx x d£ 

_d_ = % JL 

dy Cy dt 

d , d 

dt dt 


. d 

, d , d 
+ >I ' + 


( 2 . 10 ) 


In the above equations, and in those to follow, subscripts x and y, or £ and r\, denote partial differentiation 
in that coordinate direction. The only task remaining, then, is to develop expressions for the metric coef- 
ficients rj x , etc. In differential form we can write 


di = Z x dx + £ y dy + £,dt 
drj = t] x dx + rjydy + rj,dt 
di = dt 


In matrix form this becomes 
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Ut 1 U x Z y qr dx 

dt 1 = lx *ly 1, dy 

dr J [O 0 lj[rf/ 

Similarly, 

[* { x n X'ltdt 

d y = yi y n y t ^ 

dt J Lo 0 iJ[Jt 

Therefore, 

x rj x ? 

*ix n y y\ t = y ( y v y T 

_ 0 0 l J L 0 0 l 

After taking the inverse, 


’tx 

t. 

tt' 


a 


w 

* lx 

*>y 

*h 

= J 

-yf 


J'f* T - V* 

0 

0 

1 


0 

0 

iv 


where J is the Jacobian of the transformation, 

j W,n) _ f x t y 

d M *1 x >1y 

J = tx*ly ~ Zy*lx ( 2 . 11 ) 

This can be evaluated from the known physical (x, y) coordinates by noting J= \jj ~ 1 and 

j - 1 _ * { 

*«.n) y t y„ 

( 2 - 12 ) 

I he metric coefficients themselves are 

tx = Jy, 

>ix = -Jy* (2.13) 

tjy = JXt 

tt=-Xrtx~y 'rty 

t h=-x^x-y?ny 
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Unless the physical coordinates (*. j.) are defined analytically as fonchons of the computa.ronal coordinates 
the metric coefficients must be computed numencaUy. 

2.4 GOVERNING EQUATIONS IN COMPU T ATIONAL COORDINATE S 

Applying the generalized gnd transformation of the previous section to equal, on (2.1) yields 

Q t + + Qjt, + - w* ~ Ey ^ x ~ Fy & ~ ¥ ^ y = 0 (2 ' 14) 

SSgSSS&SSSSSsS 

subtracting like terms. For example, the K»sr t crm becomes 

E,<* 


M*wn 


Doing this for all the terms, and rearranging, results in 


(n 


+ 


E^x + F £y + Q St 

1 

E[/ Qx + F r £y 


t L 


F} ?x + Fr ?y "F Q*lt 


1 


j” E{/ r, x + Vy >?y 


Jx 

y 




- (F - ¥y) 


(2.15) 


The last three terms in braces, are called the metric invariant terms. By using the expressions for the metric 
strong conservation law form of the governing equations has been recovered. 


Equation (2.15) can be rewritten as 

A 

dQ 


_ + iL + iL 

dr d£ dr] 


5E„ dF v 

- + ■ 


df 


dr, 


(2.16) 


where 


A Q 


E = y (E£ x + V£y + Qt,) 
F = y (E>j x + Frj^ + Qr,t) 


3 This is not necessarily true in three dimensions, however. 
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(2.17a) 


E„ = y( E^ + F^,) 

A ] 

*v = j(£y1 x + Fytiy) 

Lsing equations ( 2 . 2 a) through ( 2 . 2 e) these can be expanded as 

Q = y 0 pu pv E T f 

P“tx + P v ty + P^t 

£ _ J_ (p“ 2 + P)ix + pUViy + put I 

PUVtx + ( pV 2 + p)ty + pvt, 

(E t + p)ut x + (E t + p)vty + E t t, 

pu^x + P Vr ty + PI t 
(pu + p)y\x + puv*ty + purl, 
puvrjx + (pv 2 + p)tf y + pvti, 

(E r + p)utj x + (E t + p)vtj y + E T »f f 

0 

r xxix + T JC yiy 
T xyi x “b T yyiy 
PxZx + Pyty 

0 

T xx r lx "b T xy*ly 
T xy r lx "b T yy t !y 

Pjtfx "b Py*ly 

where 


E„ = 


1 1 

J Re r 


A 

F>/ = 


J 1_ 

J Re, 



(2.17b) 


(2.17c) 


(2. 1 7d) 


(2-I7c) 


Px ur xx + vt X y -p— q x 

Py — Uf X y + VTyy -- -- q y 


j? ^ viscous terms, the shear stresses and heat fluxes are defined exactly as in equations ( 2 . 3 ), except 
the derivatives in the Cartesian coordinate directions must be evaluated using the chain rule. For example, 


— - J^L % . du_ 

dx di, Cjc + dr, ** 


A A A A 

Note that F and ¥ v have exactly the same form as E and E Kl but with ( replaced by > 7 . 
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3.0 TURBULENCE MODEL 


As noted briefly in Section 2,0, for turbulent flow the Reynolds stress and turbulent heat flux terms are 
modeled using the Boussinesq approach. An effective viscosity is thus defined as p = where is 

the laminar, or molecular, viscosity coefficient, and f.i ( is the turbulent viscosity coefficient. Similarly, an 
effective second coefficient of viscosity is defined as X = and an effective thermal conductivity coef- 

ficient is defined as k = k t -f k r 

The turbulent coefficients must be computed using a turbulence model appropriate for the flow being 
computed. In version 1.0 of PRO I EUS, a generalized version of the algebraic eddy viscosity model of 
Baldwin and Lomax (1978) is used to compute for wall bounded flows, (i.e., boundary layers), the 
Baldwin-Lomax turbulence model is a two-layer model, with 


(0 dinner for y n <y b 

AfiH (3.1) 

O t) outer for ^ >yb 

where y n is the normal distance from the wall, and y b is the smallest value of y n at which the values of p t from 
the inner and outer region formulas arc equal, f or free turbulent flows (i.e., mixing layers, jets, and wakes), 
Mr = (v^outer' l n inner region, in addition to the Baldwin-Lomax model, an alternate expression first 
presented by Spalding (1961), and later by Klcinstein (1967), is also available. 

In a simple boundary layer analysis, with only one solid surface, the procedure for computing fi t is rel- 
atively straightforward. In a general Navicr- Stokes analysis, however, any or all of the boundaries may be 
solid surfaces. If both boundaries in a given coordinate direction are solid surfaces, the turbulence model 
is applied separately for each surface. An averaging procedure is used to combine the resulting two f.i t 
profiles into one. If neither boundary' in a given direction is a solid surface, the formulation for free turbulent 
flows is used. In addition, values of p ( are computed separately for both the £ and r\ directions. This results 
in two complete turbulent viscosity fields. Another averaging procedure is then used to compute a single 
value of fi ( at each point in the flow. 4 

3.1 OUTER REGION MODEL 


The outer region turbulent viscosity at a given £ or rj station is computed from 

(Pi) outer ~~ F ^cpP^ Kleb^wake (3-2) 

where K is the Clauser constant, taken as 0.0168, C cp is a constant taken as 1.6, and p is the static density. 
The parameter F wake is computed from 


F wake = min < 


I ymax F max 

r v 2 — 

C wk V dlff f 


where C wk is a constant taken as 0.25, and 


(3.3) 


4 This discussion is for the most general situation. When the flow is expected to be predominantly in one direction, 
input parameters in the PROTEUS code should be used to specify that direction. 
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V diff- I I ^ I min 

where V is the total velocity vector. For wall bounded flows, | V | mi „ is taken as 0. The parameter F max in 
equation (3.3) is the maximum value of 


l\v n ) = 




for wall bounded flows 
for free turbulent flows 


(3.4) 


and y max is the value of y n corresponding to f' milx - For wall bounded flows, y n is the normal distance from 
the wall. For free turbulent flows, two values of F max and y„. x are computed - one using the location of 

| V\ max as an origin for y n , and one using the location of | V\ mm . The origin giving the smaller value of 

y max is the one finally used for computing y n , F max , and y max - In equation (3.4), | Q | is the magnitude of the 
total vorticity, defined for two-dimensional planar flow as 


Q 


dv _ du_ 
dx dy 


(3.5) 


The parameter A + is the Van Driest damping constant, taken as 26.0. The coordinate y + is defined as 


4“ 

y 


Pw“?yn 


Pw 




-y n 


(3.6) 


where u T = jrjp~ is the friction velocity, t is the shear stress, and the subscript w indicates a wall value. 
In PROTEUS, r w is set equal to ju„|£l|„. 

The function F Klei in equation (3.2) is the Klebanoff intermittency factor, given by 


I' Kieb ~ 


1 + 


< 


C Kleb yn 
ymax 



(3.7) 


where B and C KUb are constants taken as 5.5 and 0.3, respectively. This factor accounts for the exper- 
imentally observed fact that, as the free stream is approached, the fraction of time the flow is turbulent de- 
creases. 


3.2 INNER REGION MODEL 
3.2.1 Baltl» in-Lomax Model 

The inner region turbulent viscosity in the Baldwin-Lomax model is 

(dinner = P1 2 \U\ ( 3 - 8 ) 

where / is the mixing length, normally given by 

(3.9| 

and k is the Von Karman constant, taken as 0.4. 

A modified form of equation (3.9), proposed by Launder and Priddin (1973), may also be used. This 
formula is most useful for flows with steep negative gradients of shear stress normal to the wall, such as 
accelerated flows or flows with suction. I heir modified formula for I is 
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(3.10) 


where 


Pt & 

+ T ' 



and /? is a constant taken as 1.7. 

3^2.2 Spalding-KIcinstein Model 

1'he inner region turbulent viscosity in the Spalding- Kleinstein model is 

^dinner = W**~* B \* KU ~ 1 “ KU+ ~ \ ) 2 ] 


where 


V\ V 

T yf^wfPw 


(3-11) 


Again, in PROTEUS, t„ is set equal to 

3.3 AVERAGING PROCEDURES FOR MULT IPLE BOUNDARIES 

As noted earlier if both boundaries in a given coordinate direction are solid surfaces, the turbulence 
model equations are applied separately at each surface. It is assumed that the two inner re^ons do not 
overlap Vhe outer re Jons, of course, do overlap, and an averaging procedure is used to combme he two 
outer region M profiles into one. For example, if the , = 0 and r, = 1 boundaries are both solid surfaces,’ 
the two values of F wak , at a particular £ station are combmed using the following averaging formu a. 


^wake 


(Fw ake )\ /] T (Fwakellfl 

f\ +fl 


(3.12) 


Here (F w0 *,), and (F wak ,) 2 are the separate values computed for the >j = 0 and tj - 1 surfaces using equation 
(3.3). The parameters f and f 2 are defined by 


/,= 


/ 2 = 


2D, 

(y«) i 


2D 

0 


£iY 

>n ) 2 ) 


where n is a constant taken as 2.0, 0O, and (yj 2 are the normal distances to the r, - 0 and i, - 1 surfaces, 
respectively, and D x and D 2 are the normal distances from the two n surfaces to the location of | V\ max . In 
addition, the yjy m „ value needed in equation (3.7) for F KM is computed for both n surfaces and the mini- 
mum is used. These values of F WQk . and F KUb are then used in equation (3.2) to compute (/*«)«.„• 

The averaging procedure described above computes a single n, profile from the two profiles that are 
computed when both boundaries in a given coordinate direction are solid surfaces. We still must averag 


s An analogous procedure is used for solid surfaces in the { direction. 
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the twj values that result from computing /i, separately for both coordinate directions. 6 Following Goldberg 
and Chakravarthy (1987), this is done using the following formula: 8 


A*t = 


(f*t ly») i "i" ly^i 


+ ( 1 /^) 2 ] 


2 -. 1/2 


0„)2(^/)i + On)i(ju f ) 2 

[O'*)? + (yn)i ] ]l 2 


(3.13) 


Here (/*,), and (/ u,) 2 are the separate values computed due to the presence of boundaries at £ — 0 and £ = 1 
and at >7 — 0 and r, = 1, respectively. If there is only one solid surface in the £ direction, (y ), is taken as the 

toTe^lose^onf Ifth f b ?', h ^ oundari P aa ' solld surf aces, (y„), is taken as t£ normal distance 

the closest one. f there are no solid surfaces in the £ direction, (y„), is the normal distance to the location 

of either I \ I max or | as described in Section 3.1. Analogous rules are used for (j’ n ) 2 

3.4 TRANSITION MQDF.l. 

After ^ has been computed using the procedure described in the previous sections, a transition inter- 
mittcncy factor may be apphed to simulate laminar-turbulent transition. The transition model is based on 
ne given by Cebeci and Bradshaw (1984) for boundary' layer analyses, and assumes that a geometric leading 

fcr repon ,ha ' ,hc modcl is va “ d for ad,aba, ‘ c " ows at Mach 


/ft 


1 0 for x < x lr 

[ftr/ft for x > x lr 


(3.14) 


X 1S ,. thC dlst ^ nce from thc fading edge, the subscript ir indicates a value at the start of the transition 
region, and y„ is a transition intermittency factor given by 


Ytr = 1 - exp 


—G(x — x 


-■ft 


dx 


(3.15) 


In equation (3.15), u . is the velocity at the edge of the boundary layer. The factor G is given by 


G = 8.33 x 10" 


■ Re\ 


- 1.34 


where Re, u - (u t x/v)„, and v is the laminar kinematic viscosity at the edge of the boundary layer, 
written a S aSSUme ^ tHr ° Ugh the transition region, u t ~ (w e )„ and v ~ v, f , then equation (3.15) may be re- 

y f , = 1 -exp[-8.33x 10- 4 /?e^ 6 (-^- l)'] (3.16) 

Io implement equation (3.16) in PROTEUS, we replace xjx, r with ReJRe , where Re, is defined as 


Re x = 


V D 

max 


For flows predominantly in the £ direction, \ v\ mn is the maximum total velocity magnitude at the current 
S Station - D is the distance from the P oint where \v\ = \ V\ max to the leading edge at £ = 0, and v is eval- 


6 As noted earlier, this discussion is for the most general situation. When the flow is expected to be predominantly 
in one direction, input parameters in the PROTEUS code should be used to specify that direction. 
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uated at the point where | U| — | V\„, aJt - An analogous definition of Re x is used for flows predominantly in 
the y] direction. 

3.5 TURBULENT VALUES OF A AND k 

The turbulent second coefficient of viscosity is simply defined as 

) Vt ( 3 * 17 ) 


The turbulent thermal conductivity coefficient is defined using Reynolds analogy as 


k t 


c pPt 

~Pn 


(3.18) 


where c p is the specific heat at constant pressure, and Pr t is the turbulent Prandtl number. In PROTEUS, 
the turbulent Prandtl number may be treated as constant, or as a variable using the following formula 
(Wassel and Catton, 1973): 


Pr r 


Cj 


1 — exp 


Pr 3 


Q 


Pr4 


utl^t 


C Pr\ Pr l 


c. 


“n-7^7) 


(3.19) 


Here C PrU C Pr2 , C PrZ) and C Pr4 are constants taken as 0.21, 5.25, 0.20, and 5.0, respectively, and Pr { — 
is the laminar Prandtl number. 
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4.0 TIME DIFFERENCING 


The governing equations are solved by marching in time from some known set of initial conditions using 
a finite difference technique. The time differencing scheme currently used in PROTEUS is the generalized 
scheme of Beam and Wanning ( 1978). The time derivative term in equation (2.16) is wntten as 


A A 

dQ AQ” 

dr At 


0, 3(AQ") 

1+02 


+ 0 


(•' 


J_ 

2 



02 

1 + 0 2 


AQ 


n— 1 


At 


or, 


AQ" = 


0, At 

1 + 0 , 


d(AQ") 

dz 


+ 


At 

1+0, dr 


+ 


1 + 0 - 


AQ 


n— 1 


o[(o, - ~ - 0 2 ) ( At) 2 + (At) 3 ] 


(4.1) 


where AQ n = (2 ' 1 — Q". The superscripts n and n + 1 denote the known and unknown time levels, re- 
spectively. 

The parameters 0, and 0 2 determine the type of time differencing scheme used. Some of the methods 
available with the above formula are given in the following table. 


0. 

0 2 

Method 

Truncation Error 

0 

0 

Euler explicit 

6>(At) 2 

0 

-1/2 

Leapfrog explicit 

O(At) 3 

1 

0 

Euler implicit 

O(At) 2 

1/2 

0 

Trapezoidal implicit 

0( At) 3 

1 

1/2 

3-point backward implicit 

O(At) 3 


Note that even though the generalized time differencing formula includes explicit methods, the I RO 1 LUS 
code assumes an implicit method is being used. Note also that the truncation error listed in the table is the 

error in the expression for AQ '. The overall numerical method used in modelling the differential equations 

requires AQ"/At, so the order of the overall method is this truncation error divided by At. 

Solving equation (2.16) for dQ/dr and substituting the result into equation (4.1) for d(AQ")/dx and 
dQ"/dr yields 


A„ 0! At ( d(AE”) d(AF”) \ At ( dE” , d¥* \ 

AQ = -7T0^^^ + ^r J-t T 0r(,-dr + ^rj 

0, At / d(AE^) d(AE^) \ Az f dE^ dl> \ 

+ l+0 2 l dl + dn J + 1 +0 2 ^ dl dr, J 

+ AQ"" 1 + 0^0, - y - 0 2 )(At) 2 + (At) 3 J (4.2) 
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5.0 LINEARIZATION PROCEDURE 


5.1 INVISCID TERMS 


Equation (4.2) is nonlinear, since, for example, AE* = E" +l — E" and the unknown E 7 ^ 1 is a nonlinear 
function of the dependent variables and of the metric coefficients resulting from the generalized grid trans- 
formation. The equations must therefore be linearized to be solved by the finite difference procedure used 
in PROTEUS. This is done by expanding each nonlinear expression in a Taylor series in time about the 
known time level n. Letting G represent any nonlinear expression, 

G n+ ' = G n + At + 0(A?) 2 (5.1) 

where 

dG = dG dp qg d(pu) dG d(pv) dG dE r 

St dp dr + d(pu ) dr + d(pv) dr + dE T d-z 


Note that for linearization purposes only the metric scale coefficients have been assumed to be locally inde- 
pendent of time. Note also that for this linearization procedure to be second-order accurate, dG/dr (and 
therefore dpjdr , d(pu)ldr y etc.) need only be first-order accurate. Using forward differences, then, so that 



k+i 


At 


• + O(At) 


a n 

At 


T O(At) 


etc., equation (5.1) becomes 

c ” +i - «■ + ( t ) v + ( & y** + ( m y^ + 


(5.2) 


As an example the d(puv£ y )jd< l term from the ^-momentum equation (part of the second element of 

dE/<3£) will be used. The nonlinear part of this term is (puv) n+l . Rewriting this in terms of the dependent 
variables, 


(puvf+ ] 


( P u)(pv ) t +1 
P 


Using equation (5.2), this is linearized as 

(pwv) n+I = (puv) n - (w)”(p" +l - p n ) + v”[(pw)" +1 - fcu)"] + «”[(pv)" +1 - (pv)”] + 0{ At ) 2 
which can be rewritten as 


A(puv) n — — (uv) n Ap n -f v n A(pu) n -f u n A(pv) n + O(At) 2 


This linearization procedure, when applied to the entire AE n term in the vector equation (4.2), can be 
written as 
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where (dEjdQ) n is a Jacobian coefficient matrix (not to be confused with the Jacobian J of the generalized 
grid transformation.) A similar equation can be written for AF n . 

A A 

Each term in each element of E and F, given by equations (2.17b) and (2.17c), is linearized using the 

A A A A 

above procedure to generate the elements of the Jacobian coefficient matrices dEjdQ and dEjdQ. (Note 
that dEjdQ = JdEjdQ.) When this is done dEjdQ can be written as 


dE 

A 

3Q 



ix Zy 0 

dp dp dp 

dp „ dp dp 

f2 ^ x+fi ~Su) 


(5.4) 


where f = u£ x + v£ y and f 2 = (E r + p)jp. The Jacobian matrix dEjdQ has the same form as dEjdQ , but 
with £ replaced by rj. 

The linearized pressure terms have deliberately been left in terms of dpi dp, dpjd(pu), etc. The ex- 
pressions to be used for these derivatives depend on the equation of state. Those currently built into the 
PROTEUS code, for a perfect gas, are presented in Section 5.3. 

5.2 VISCOUS TERMS 


The nonlinear viscous terms in equation (4.2), involving AE n K and AF'J,, must also be linearized. To do 

A A 

this, the elements of E v and E V} given in equations (2.17d) and (2.17e), must first be rewritten in terms of 
the dependent variables, and with derivatives in the Cartesian directions transformed to derivatives in the 
computational directions using the chain rule. When the resulting expressions arc substituted into equation 
(4.2), mixed second derivatives appear as well as second derivatives in a single coordinate direction. The 
mixed, or cross, derivative terms would lead to considerable complications in the implicit numerical solution 
algorithm if they were linearized using the procedure presented in Section 5.1. The two types of second 

A A 

derivatives are thus treated differently, and E v and F^ are written as 


AAA 

Ej/ = E^ -I - 


F y = F + F 


Vl 


(5.5) 


A A ^ A 

where E v and F y only contain derivatives in the £ and >/ directions, respectively, and E y 2 and E Vl contain 

derivatives in the other direction. The fully expanded expressions for E K[ , E Vr etc., are fairly long, and 
therefore are presented in Appendix A. 
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5.2.1 Non* Cross Derivatives 


Examination of the elements of E V] in equations (A. 2a) through (A.2c), and (A.2e), shows that every 
term has the form fg ( , where g is a function of the dependent variables, and/is a function of p, k, k, and/or 
the metric coefficients. Expanding in a Taylor series about time level n gives 


(/g / +1 



d{fg>) - 

di 


n 

At + 0( At) 2 


For linearization purposes only, we will assume f is locally independent of time. We can thus write 

(/^r 1 = (/gsf +r -§z [ % ] ^ + o(at ) 2 


where 


Therefore 


dg _ dg dp^ dg d(pu) + 

dr dp dr d(pu) dr 


(fg t ) H+l 


= (/g t ) n +r 


d_ 

di 


dg_ 

dp 


A p + 


dg 

d{pu) 


A (pu) + - 


+ O(At) 2 


As with the inviscid terms, the linearization procedure for the entire AE^ viscous term in equation (4.2) can 
be written as 


AE 


n 

V\ 


dE L 


dQ 


AQ” + 0(At) 2 


A similar equation may be written for AFfy 1 he Jacobian coefficient matrix d¥. y JdQ is 


(5.6) 



PROTEUS 2-D Analysis Description 


Linearization 27 



a yy ~ + ( 2/1 + 2 )^ 

a xy = (/‘ + 2)^ 



Like the pressure terms discussed earlier, the form of the temperature terms will depend on the equation 
of state being used. Those currently built into the PROTEUS code, for a perfect gas, are presented in 
Section 5.3. 


Note that in equation (5.6) the derivatives appearing in the Jacobian coefficient matrix dE v jdQ are also 

A l 

to be applied to the AQ 71 appearing outside the parentheses. For example, the element in the second row 

A A 

and second column of dE y JdQ , which corresponds to the A (pu) term in the r-momentum equation, is 
a XJt d(l/p)/df . For this term, the notation used in equation (5.6) means 

n d ( \ 

= « 


The Jacobian coefficient matrix for the remaining non-cross derivative viscous terms, 5F K /3Q, has the 

A A 1 

same form as dEyJdQ, but with f replaced by r\. 

5.2.2 Cross Derivatives 

As stated earlier, linearizing the cross derivative viscous terms in the same way as the remaining terms 
is very complicated within the framework of the implicit numerical solution algorithm used in PROTEUS. 
They are therefore simply lagged (i.e., evaluated at the known time level n and treated as source terms.) 
As noted by Beam and Warming (1978), this does not lead to a formal accuracy loss since 
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(5.8) 


AEy 2 = AEy~' + Q(A t ) 2 
AF*. = AF^T 1 + 0 { At) 2 


5.3 EQUATION OF STATE 


The expressions to be used for dpjdp , dTjdp , etc., which arise from the linearization procedure, depend 
on the equation of state. The equation currently built into PROTEUS is for perfect gases, and can be 
written as 

p = (y- l)[£ r --±-p(i/ 2 + v 2 )] (5.9) 

or, in terms of temperature, as 

7 '-i[ J ! : -T (“ 2 + v2 )] < 51 °> 


With this equation of state, then, the appropriate derivatives are 


dp_ 

dp 



d_P_ 

d(pu) 


= - (y - i)w 


d{pv) 


= - (y - l)v 


dr 

dp 


d _P _ . 

dE T y 



dT u 

d{pu) C \P 

dT v 

d(pv) c vP 

dT _ 1 
dEj c vP 


(5.11a) 
(5.11b) 
(5.11c) 
(5. lid) 
(5.12a) 

(5.12b) 
(5.12c) 
(5.1 2d) 


If constant stagnation enthalpy is assumed, as discussed in Section 2.2, the appropriate equation of state 
is 

P = ^~Y~p\^T-\ (“ 2 + v2 )] ( 513 ) 

and the temperature becomes 

t =W t ~ t ( “ 2+v2) ] (514) 

With these equations, the derivatives of p and T with respect to the dependent variables are 
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(5.15a) 


dp 


i t ^[ / v+t ( w2 + v2 )] 

dp y - 1 

d{ P u) y 

dp y - 1 

a(pv) y 


dT 

dp 


c pP 


1 (w 2 + *' 2 ) 


ar _ u 

d(pu) c pP 


dT _ v 
d(pv) c p/> 


5.4 LINEARIZED GOVERNING EQUATION 

The linearized form of equation (4.2) can now be written as 



n— 1 


1 + 0 


2 AQ" -1 + O 


(0,-y-0 2 ) 


(At) 2 , (0 3 — 0|)(At) 2 , (At) 3 


(5.15b) 

(5.15c) 

(5.16a) 

(5.16b) 

(5.16c) 


(5.17) 


There are a couple of things that should be mentioned about this equation, first, this equation is in 
so-called "delta" form. We will actually be solving this equation for A Q n and recovering Q ,+l from 
q«+i = A O' 1 + Q" . And second, in the coefficients of the cross derivative viscous terms the time differencing 
parameter 0, has been replaced by 0 3 . For second-order time differencing (i.c., if 0, = 0 2 + 1/2), 0 i should 
be set equal to 0,. For fust-order time differencing, however, 0 3 can be set equal to zero without losing 
accuracy. 
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6.0 SPACE DIFFERENCING 


To solve equation (5.17) an evenly spaced grid is defined in the computational ({, rj) coordinate system. 
Spatial derivatives are then approximated by finite difference formulas. First derivatives in the £ direction 
are approximated using the following variably centered formula. 



I — &(flj “ *TT [(^ — a )fi+\J + ( 2a \)fij ^ 

U ^ 


( 6 . 1 ) 


where a = 0 for first -order forward differencing, a = 1/2 for second-order central differencing, and a — 1 for 
first-order backward differencing. All of the PROTEUS cases run to date have used central differencing. 
The subscripts / and j represent grid point indices in the £ and rj directions. The computational grid spacing 
is constant, and equal to 1 /(jV| — l), where is the number of grid points in the £ direction. A similar 
formula is used for first derivatives in the yj direction. 


The non-cross derivative viscous terms in the £ direction in equation (5.17) all have the form 


where Q represents one of the elements of Q. Using central differences this is approximated by 


~k l/w fe4 ®] ( “ 

(Af) 

-ft- i/2.y[(?A0 U -(gA0 t _ 1>y ]} 

= — ^ {(4; + y; + i,y)[0? A C) i+ .,y - 

2(A0 

- (4y + ft- 

= l -T «4-i.y + 

2(A 0 

— ( f'l— 1 j y + 2fi J 4- j)C?^0/,y 

+ (4y+ W^W ( 6 - 2 ) 

A similar formula is used for second derivatives in the >7 direction. 

Cross derivative viscous terms are evaluated using the following central difference formula. 


PROTEUS 2-D Analysis Description 


Space Differencing 31 



d l 


, ^ hi 


r- (. ?): ; 


1 


^ ^ * 1 - /■ ^7 H J -h - 1 .yl 1 ,y] 

dA£A^~ ^‘ + h/&+ 1 J+i ~~ 1 .y-i ) 

” Ji--\ji%i-\j+] ~~ Si- I ; y_ 1 )] 


(6.3) 


Note that this formula is only needed for the source terms, since the viscous cross derivative terms are 
lagged. 


When first derivatives are needed normal to a computational boundary, such as for Neumann boundary' 
conditions, either first- or second-oider one-sided differencing is used. The first-order formula at the £ = 0 
boundary' is 


Of \ i 

dt J, J~ aT ( ^2j- J\j) 


and at the Z — 1 boundary, 


df 


d£ Js j~ A£ ^ v i -j ^1-1 .y) 
The second-order formula at the Z = 0 boundary is 

2S{-<-Vw' + 4 /2J-/3,y> 

and at the £ = 1 boundary, 


J 

dZ 


1 J 


oz 


AW 


2A£ 


( /v, - 2 ,j 4 /v, - 1 J + 3/v, j) 


(6.4) 


(6.5) 


(6.6) 


(6.7) 


Similar formulas are used at the y - 0 and >; — 1 bomularics. 
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7.0 BOUNDARY CONDITIONS 


Choosing boundary conditions is perhaps the most important step in solving a flow problem with 
PROTEUS. Since the equations being solved at interior points are the same for every problem, the 
boundary conditions are what determines the final flow field for steady flows. 

With the difference formulas presented in Section 6.0, N eq boundary conditions are required at each 
computational boundary, where N eq is the number of equations being solved. Note, however, that this is 
a numerical requirement, not a mathematical one. For example, for one-dimensional Euler flow N tq = 3. 
However, characteristic theory shows that, mathematically, only two conditions may be specified at a sub- 
sonic inliow b- r.;i:d::r /, and only one at a subsonic outflow boundary (Pulliam, 1986 a). Some sort of ex- 
trapolation is typically used for the additional numerical boundary conditions. 

A variety of boundary conditions are built into the PROTEUS code, including: (1) specified values 
and/or gradients of Cartesian velocities u and v, normal and tangential velocities V n and V n pressure p , 
temperature T, and density p\ (2) specified values of total pressure p T , total temperature T t , and flow angle; 
(3) linear extrapolation; and (4) spatial periodicity. Another useful boundary condition is a "no change from 
initial condition" option for u, v, p, T y p, p T) and/or T r . Provision is also made for user- written boundary 
conditions. The boundary conditions may be steady, unsteady, or time-periodic. The exact combination 
of boundary conditions to be used will depend on the problem being run. 

The boundary conditions in PROTEUS are treated implicitly. They may be viewed simply as additional 
equations to be solved by the ADI solution algorithm. And, in general, they involve nonlinear functions 
of the dependent variables. They must therefore be linearized using the procedure described in Section 5.0. 
The following sections describe this linearization for the general types of boundary conditions currently built 
into PROTEUS. 

7.1 NO CHANGE FROM INITIAL CONDITIONS, Ag = 0 

This boundary condition simply sets the boundary value of the function g equal to its initial condition 
value. It can be written as 


a n fl+1 n n 

=g -g = o 


(7.1) 


In general, g can be a nonlinear combination of the dependent variables Q. Linearizing g using the proce- 
dure described in Section 5.0, we get 


n+1 n , / d g 

g *=*+( — 

V dQ 


AQ" + O(Ar) 2 


(7.2) 


Neglecting the 0 (At) 2 Linearization error, the linearized form of equation (7.1) can thus be \yritten as 

( dg ^ " 


A 

dQ 


aq ,! = o 


(7.3) 


7.2 SPECIFIED FUNCTION, 

A specified function at a boundary can be written simply as 


g n+X =/ 


(7.4) 
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where g is the function being specified and /is the value being specified. Note that /can vary along the 
boundary, and can be time-dependent. Using equation (7.2) and neglecting the linearization error, the 
linearized boundary condition becomes 

^JaQ"=/-/ (7.5) 


7.3 SPECIFIED COORDINATE DIRECTION GRADIENT, d.e/d</> = [ 


A specified gradient of a function in a coordinate direction can be written as 



= / 


(7.6) 


where g is the function whose gradient is being specified, /is the specified value, and <f> is the coordinate 
direction f or rj. Note that / can vary along the boundary, and can be time-dependent. 

The linearized form of g is given by equation (7.2). The linearized form of equation (7.6) can thus be 
written as 




A 

dQ 



= /+ O(At) 2 


(7.7) 


Replacing differential operators with difference operators and neglecting the linearization error, the 
linearized boundary condition can be written as 



=/- v” 


(7.8) 


where </ represents the one-sided difference operator to be used at the boundary. Options are available in 
PROTEUS to use either first-order two-point or second-order three-point differencing. 


Note that this boundary condition is a specified value of the derivative with respect to the computational 
coordinate, not with respect to the physical distance in the direction of the computational coordinate. 
Following Korn and Korn (1968), and using the properties of the generalized coordinate transformation, it 
can be shown that for the £ direction the two derivatives are related by 


dg 


ds 


£ 


J 


V 


/ n\ + 


dg_ 


Similarly, for the tj direction, 


Sg _ 7 dg 

N VdT? d ” 


If the value /= 0, of course, the two derivatives are equivalent. 

7.4 SPECIFIED NORMAL DIRECTION GRADIENT, Vg.n =/ 

A specified gradient of a function normal to the boundary can be written as 

V +, -«=/ (7.9) 

where g is the function whose gradient is being specified, /is the specified value, and n represents the unit 
vector normal to the boundary. Note that / can vary along the boundary, and can be time-dependent. 
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For illustrative purposes, assume wc are specifying a gradient normal to a constant f boundary. Then 


n = - 


|V{| 


■ = -*-£ 
m 


/ + 






where 

" 2 = Vd + ^ 


Equation (7.9) can then be written as 

-^fe£ +1 {x + *y +1 £y)=/ 

Using the chain rule to expand g^ x and 



Substituting into equation (7.10) and rearranging, 

g{ + '(£ + (}) + gq + \t*x + ty*y) = m f 


(7.10) 


Solving for gV x , 



m 


T £y*1y) 


m 


dg_ 

drj 



(7.11) 


Now, in order to incorporate this equation into the ADI solution procedure used in PROTEUS, the 
dg/drj term in equation (7.11) is lagged one level, and evaluated at time level n instead of n + 1. Strictly 
speaking, this introduces an O(At) error into the solution. In practice, however, the actual error will depend 
on the degree of nonorthogonality of the coordinates near the boundary. For orthogonal coordinates no 
error is introduced. 


Using equation (7.2), and introducing difference operators and neglecting the linearization error, we can 
now write the linearized boundary condition as 


St 





+ W V” “ 6 S g ” 

m 


(7.12a) 


where 3 { represents the one-sided difference operator to be used at the boundary. Options are available in 
PROTEUS to use either first-order two-point or second-order three-point differencing. 


Specifying a gradient normal to a constant rj boundary is done in an exactly analogous manner. The 
resulting equation is 
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> / 



A 

dQ 



/ 

m 


m 


(7.12b) 


where 


r 2 . 2 
m = v Vx + Vy 


7.5 LINEAR EXTRAPOLATION 


linear extrapolation from the two adjacent interior points is also available as a boundary condition. 
At the £ = 0 boundary 7 , where 1=1, this can be written as 
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Analogous extrapolation boundary conditions can easily be written for the remaining boundaries. 
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8.0 SOLUTION PROCEDURE 


8.1 ADI ALGORITHM 

I hc governing equations, presented in linearized matrix form as equation (5.17), are solved by an alter- 
nating direction implicit (ADI) method. 1 he lorm of the ADI splitting is the same as used by Briley and 
McDonald (1977), and by Beam and Warming (1978). Although the split equations can be developed in 
more than one way, in this discussion the approximate factorization approach is used. 


Letting I 1 IS(5. 1 7) represent the left hand side of equation (5.17), we can write 


LHS(5.17) 




( 8 . 1 ) 


where I represents the identity matrix. Note that in this equation, using the djd £ term as an example, the 
notation used is meant to imply 


a 


dE dEy i 


dQ dQ 


AQ — ~rr [ #AQ-- ^-AQ 


A 

dE 


v. a 


dt, 


dQ 


3Q 


The term in braces in equation (8.1) can be factored to give 


LHS(5.17) = 


1 + 


At g / dE 


d Ei 


i + a<j 


-( 
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dQ dQ 
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dQ 
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ar 

A 

dQ 


dQ 


AQ 


(8.2) 


The last term represents the splitting error. Note that, since AQ" = O(Ar), this term can be neglected 
without affecting the overall time accuracy of the algorithm, even when second-order time differencing is 
used. 
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Equation (5.17) can thus be rewritten in spatially factored form, and, neglecting the temporal truncation 
and splitting error terms, becomes 


1 + 


0i Ar 


dE 


dE. 
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dQ dQ 
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aq” = 


V (1+®3)At ( <& v 
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V 2 
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dr, 


(8.3) 


Equation (8.3) can be split into the following two-sweep sequence. 
Sweep 1 (£ direction) 


AQ + 


g i Ar d 

l+0 2 df 


a \ n 

dE 
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AQ 


°i At d 
1 + 0 2 d{ 
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+ (1 +0 3 ) A t ( ^V 2 Y 0 3 At 
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dZ drj 
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dE v 
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1 + 0 - 


AQ” _1 (8.4a) 


Sweep 2 (u direction) 


AQ" + 


°i At 0 
1 + 0 2 dr, 


dF 

A 

0Q 


AQ 


#i At d 

1 + 0 2 dr, 


dF. 


dQ 


AQ” 


A Q 


(8.4b) 


In the above equations, Q represents an intermediate solution to the governing equations. 7 It should be 
noted that in PROTEUS, physical (i.e, n + 1 level) boundary conditions are used during the first ADI 
sweep. This introduces an O(At) error in dQjdr on the boundary for unsteady flows, but no error for steady 
flows. This point is discussed in detail by Briley and McDonald (1980). 


The notation here is somewhat inconsistent. The quantity AQ" = Q"*' — Q", but AQ' = Q' — Q", not Q n+1 - Q'. 
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Applying the spatial differencing formulas of Section 6.0 results in 


Sweep 1 (£ direction) 


AQ, + 


At 
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(8.5a) 


Sweep 2 (*; direction) 


a 0 j At 
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aq 


(8.5b) 


The subscripts i and y represent grid point indices in the £ and rj directions. For notational convenience, 
terms without an explicitly written / or j subscript are understood to be at i or y\ In the viscous terms on 

the left hand side, / is the coefficient of d/d£ (or d/d*/, depending on the sweep) in the dE Ki /dQ (or 
d¥ yJdQ) Jacobian coefficient matrix. Similarly, g is the term in the parentheses following d/d£ (or d/d*/) 

in the dE Ki /dQ (or dF Kj /dQ) Jacobian coefficient matrix. Equations (8.5a) and (8.5b) represent the two- 
sweep alternating direction implicit (ADI) algorithm used to advance the solution from time level 
n to n 4- 1 . 


8.2 MATRIX INVERSION PROCEDURE 
8.2.1 Non-Periodic Boundary Conditions 


The complete set of algebraic equations for the first ADI sweep with non-periodic boundary conditions 
can be written in the following block matrix form. 8 


8 Although this discussion is written for the first ADI sweep, an exactly analogous procedure is followed for the 
second sweep. 
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- 

A * 


- 

B', C'l A', 

AQ, 


S', 


A* 



a 2 b 2 c 2 

aq 2 


^2 


A * 



A3 B3 c 3 

aq 3 
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S3 

• 

■ 

m 

= 

• 

■ 

■ 

A * 


• 

A,v, -2 ®/v, -2 C/v, _ 2 

AQ*, -2 


S*, -2 

A,V, -1 ®AT, -1 C„,_, 

AQ*,-, 


s*, -, 


A * 



C'v, A' V| B'*, 

AQ*, 


S 'v, 

- 



- 


These equations result from the application of equation (8.5a) for / = 2 to N x — 1, with boundary conditions 

added at / = 1 and i — N v The parameter AQ’ is the element vector containing the unknown dependent 
variables; A, B, and C are the N tq xN ttl coefficient submatriccs at / — 1, /, and i + 1, respectively; and S is the 
jV^-element subvector containing the explicit source terms. Also, A', B', and C' are the coefficient sub- 
matrices and S' the source term subvector for the boundary conditions. A variety of boundary conditions 
may be used. They are described briefly in Section 7.0, and in greater detail in Volumes 2 and 3. 


Note that the equations at the boundaries may contain coefficients at the boundary point and the two 
adjacent interior points. This occurs, for example, when extrapolation or second-order gradient boundary 
conditions are specified. As written, therefore, the coefficient matrix in equation (8.6) is not block 
tridiagonal. However, A\ can be eliminated by multiplying the second row of the matrix by A'i Cj 1 and 
subtracting from the first row. C'^ can be eliminated in a similar manner. Doing this, we define 

B, = B, - A', CJ’Aj 


c, = C', — A', c 2 b 2 

( 8 . 7 ) 

S, = S', - A', 


A*, = A'/y ( - C' v , A /V| ' _,B jV( 


B*, = BV, - C'v, a; v ;_,c^_, 

(8.8) 

S*, = S'v, — C'v, A ( v, 1 _,S*, 
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The set of algebraic equations solved during the first ADI sweep can now be written as 
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Since the coefficient matrix is now block tridiagonal, the equations can be solved using the block matrix 
version of the Thomas algorithm (e.g., see Anderson, Tannehill, and Pletcher, 1984). The procedure can 
be summarized as follows: 

1. Define D t = B r 

2. Compute E l = Df 1 C l and AQJ = Dj^Sj. 

3. For / = 2 to N { , compute 
D; = B, - 

E^D-'C, 

aq; = DT\s, - a,aq;_, ) 

(Actually, E, is only needed for / = 2 to N { — 1). 

4. Then, set AQ^ = AQ^ . 

5. Finally, for / = — 1 to 1, compute AQ, = AQJ — E,AQ 1+l . 

In the PROTEUS code, in step 2 Ej and AQ{ are actually obtained by solving = Q and 

A A 

d.aq; =s, using LU decomposition of D. A similar procedure is used to compute E, and AQ' in step 
3. 

8.2.2 Spatially Periodic Boundary Conditions 

In computational coordinates a spatially periodic boundary condition in the £ direction may be repres- 
ented as shown in Figure 8.1. 9 


9 As in Section 8.2.1, this discussion is written for the first ADI sweep, but an exactly analogous procedure is followed 
for spatially periodic boundary conditions in the second sweep. 


PROTEUS 2-D Analysis Description 


Solution Procedure 41 
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V 



Figure 8.1 - Spatially periodic boundary condition. 


The grid points along the / = 1 and i = ,V, lines are "similar" in the geometric sense, and have the same 
flow solution. Therefore, for a spatially periodic boundary condition in the £ direction, Q t — Q.Vj ■ 

To implement this boundary condition, an additional set of points is added at i=N l + 1, setting 
* y +1 = Q 2 . This allows us to use central differencing in the £ direction at i = N u computing the coefficients 
the same way as at the interior points. 


i 
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The resulting set of algebraic equations will consist of N { - 1 equations (for /= 2 to N v ) t with N t + 1 
unknowns. The block coefficient matrix thus has 1 rows and + 1 columns, as follows: 
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These equations result from the application of equation (8.5a) for /— 2 to A/^. As in the previous section, 

the parameter AQ* is the ^-element vector containing the unknown dependent variables; A, B, and C are 
the N eq xN eq coefficient submatrices at /— 1, i f and / + 1, respectively; and S is the A^-element subvector 
containing the explicit source terms. 

Since Q, = and Q 2 = 4lt equation (8.10) can be rewritten with N { - 1 unknowns as: 
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An efficient algorithm to solve this system can be derived that is similar to the Thomas algorithm for 
block tridiagonal systems. The procedure can be summarized as follows: 

1. Define D 2 = B 2 and F 2 = C^. 

2. Compute E 2 = D 2 »C 2 , G 2 - D 2 ] A 2 , and AQ 2 = D 2 S 2 


PROTEUS 2-D Analysis Description 


Solution Procedure 43 


3. For / = 3 to A 7 , — 1, compute 
D< = B, — A ( K ( _, 

E^D-'C, 

F< = ~ E/_j 
G; = - Dr' A^, 

AQ' =D- , (S / -A / AQ;_ 1 ) 

4. Compute 

G,v, -i = -t(G ; v, -i ~ A/e, -]G, Vl -2) 

*“ A', -1 = A,v, — F jV] _ 2 K v . _ 2 

/V,-l 

D v, = B ;Vl - X F A 

i=2 

F/Aq; j 

5. Then, set AQ Aj = AQ' Vi . 

6. Compute AQ Vi - AQ Aj - G Aj _,AQ^. 

7. Finally, for i =N X - 2 to 2, compute AQ, = AQ' - E ( AQ iM - G,AQ Ar 

In the PROTEUS code, in step 2 E 2 , G 2 , and AQ 2 are actually obtained by solving D 2 E 2 = C 2 , 
D 2 G 2 = A 2 , and D 2 AQ 2 = S 2 using LU decomposition of D. A similar procedure is used to compute E ( , 
G ( , and AQ' in step 3, and G Ai _j and AQ Al in step 4. 

8.3 UPDATING BOUNDARY VALUES 

8.3.1 Non-Periodic Boundary Conditions 

With the ADI algorithm described in Section 8.1, if gradient or extrapolation boundary conditions are 
used for the first sweep, the boundary values from the first sweep must be updated after the second sweep. 
This point is easiest to illustrate by looking at the following figure. 
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Figure 8.2 - Updating boundary values for non-periodic 
boundary conditions. 

In Figure 8 2 a 5x5 grid is shown in computational space. The triangles represent grid points at which 
the intermediate ’values Q' are computed during the fu st ADI sweep. These include the boundary points 
nt f. _ o an A i - i The circles represent grid points at which the final values Q" +1 are computed during ; t & 

i 

final values at the interior points. 

To do this, after the second sweep the boundary condition equations are rewritten and solved at the { 
boundaries. At the £ = 0 boundary, 

IV"AQ" + C'/’AQ" + AfAQ” = S'," ( s - 1 2 ) 

A 

All the terms in equation (8.12) are known except AQf. Solving, 

AQ, = (B',")" 1 (S'," - C'"AQ" - A',"AQ^) ( 8 - 1 3 ) 


At the £ = 1 boundary', 


C'y, AQ" W| _ 2 + Aft, AQft, + Bft, AQ" V , = Sft 
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Finally note from Figure 8.2 that new comer point values are never computed in the solution algorithm. 
To make the corner values consistent with the rest of the flow field, in PROTEUS the comer values of 

fnh ih '°! C r Crgy Ll arc arbitrarily defmed by linearly extrapolating from the two adjacent points 

in both the £ and t, directions, and averaging the two results. The comer values of the velocities are updated 
> oing the same type of extrapolation. Instead of averaging, however, the extrapolated velocity whose 
absolute value is lower is used. This was done to maintain no-slip conditions at duct inlets and exits. 

8.3.2 Spatially Periodic Boundary Conditions 

bouX'condtont^ uS" ^ ' he ^ ^ * compli “ ted »*“» periodic 
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I igurc 8.3 - Updating boundary values for periodic 
boundary conditions in the £ direction onlv. 


The situation for a periodic boundary condition in the £ direction but not in the r, direction is shown 
in Figure 8.3. I he triangles again represent grid points at which intermediate values are computed and the 
circles represent gnd pomts at which final values are computed. As can be seen from the figure the inter- 
mediate values at £ - 0 must be updated after the second sweep to be consistent with the final values at the 
interior points. This is easily done by setting Qj = for j = 1 to N 2 . 
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Figure 8.4 - Updating boundary values for periodic 
boundary conditions in the tj direction only. 


The situation for a periodic boundary condition in the r\ direction but not in the £ direction is shown 
in Figure 8.4. In this case, the intermediate values at £ = 0 and at ^ = 1 must be updated after the second 
sweep. To do this, the same procedure described in Section 8.3.1 for non-periodic boundary conditions is 

used, but for j=2 to N 2 instead of N a - 1. Then, for the lower comer values, Q,., = Qi.iv, ^d 

A A 
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Figure 8.5 - Updating boundary values for periodic 
boundary conditions in both the £ and rj directions. 


And finally, the situation for periodic boundary conditions in both the £ and rj directions is shown in 
Figure 8.5. Like the case with periodic boundary conditions only in the £ direction, the intermediate values 

at £ = 0 must be updated after the second sweep. This is again done by setting Q, = Q V] for j = 1 to N 2 . 
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9.0 ARTIFICIAL VISCOSITY 


With the numerical algorithm of Section 8.0, high frequency nonlinear instabilities can appear as the 
solution develops. For example, in high Reynolds number flows oscillations can result from the odd-even 
decoupling inherent in the use of second-order central differencing for the inviscid terms. In addition, 
physical phenomena such as shock waves can cause instabilities when they are captured by the finite dif- 
ference algorithm. Artificial viscosity, or smoothing, is normally added to the solution algorithm to suppress 
these high frequency instabilities. Two artificial viscosity models are currently available in the PROTEUS 
computer code - a constant coefficient model used by Steger (1978), and the nonlinear coefficient model 
of Jameson, Schmidt, and Turkel (1981). The implementation of these models in generalized 
nonorthogonal coordinates is described by Pulliam (1986b). 

9.1 CONSTANT COEFFICIENT ARTIFICIAL VISCOSITY 


The constant coefficient model uses a combination of explicit and implicit artificial viscosity. The 
standard explicit smoothing uses fourth-order differences, and damps the high frequency nonlinear insta- 
bilities. Second-order explicit smoothing, while not used by Steger or Pulliam, is also available in 
PROTEUS. It provides more smoothing than the fourth-order smoothing but introduces a larger error, 
and is therefore not used as often. The implicit smoothing is second order and is intended to extend the 
linear stability bound of the fourth-order explicit smoothing. 

The explicit artificial viscosity is implemented in the numerical algorithm by adding the following terms 
right hand side of equation (8.5a) (i.e., the source term for the first ADI sweep.) 

£ (2) AT £ (4 ) At 

(V^,Q + V^Q) - -£j— [(V^) 2 Q + 0VV 2 Q] (9.1) 

and are the second- and fourth-order explicit artificial viscosity coefficients. The symbols V 
are backward and forward first difference operators. Thus, 

W = OrQi-i 
* t Qi = Q M -Qi 
V ? A ? Q i = Q, +1 -2Q^Q i _ 1 

(VjAj) 2 Q / = Q, +2 - 4Q /+1 + 6Q, - 4Q / _ 1 + Q /_ 2 
Equivalent formulas arc used for differences in the rj direction. 

A few details should be noted at this point. First, the sign in front of the artificial viscosity term being 
added to equation (8.5a) depends on the sign of the "i" term in the difference formula. For damping, that 
term must be negative when added to the right hand side of the equations (i.e., explicit artificial viscosity), 
and positive when added to the left hand side (i.e., implicit artificial viscosity.) See Anderson, Tannehill, 
and Pletcher (1984) for details. Second, the terms being added are differences only, and not finite difference 
approximations to derivatives. They are therefore not divided by A£, etc. Third, the variables being dif- 

A 

ferenced are Q, not Q. As noted by Pulliam (1986b), scaling the artificial viscosity terms by 1 jJ makes them 
consistent with the form of the remaining terms in the equations. Fourth, the terms are also scaled by At. 
This makes the steady state solution independent of the time step size (Pulliam, 1986b). And finally, note 
that the fourth-order difference formula cannot be used at grid points adjacent to boundaries. At these 
points, therefore, the appropriate fourth-order term in expression (9.1) is replaced by a second-order term. 
Thus, for points adjacent to the £ = 0 and £ = 1 boundaries, — e^ > Ar[(V 4 A^) 2 Q]/y is replaced by 


to the 


where 
and A 
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(9.2) 
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4^ 

j 




A similar expression is used at points adjacent to the t, = 0 and t, = 1 boundaries. 

The implicit artificial viscosity is implemented by adding the following terms to the left hand side of the 
equations specified. 


E/At a, 

— — [V^Aj(7AQ )] to equation (8.5a) 


E/ar a 

— ~j— [^A^v/AQ )] 


to equation (8.5b) 


(9.3) 


Note that the addition of the artificial viscosity terms, in effect, changes the original governing partial 
differential equations. At steady state, the difference equations with the artificial viscosity terms added ac- 
tually correspond to the following differential equations. 10 
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The implicit terms do not appear, since they difference AQ, and in the steady form of the equations 

AQ = 0 . The artificial viscosity terms do not represent anything physical. The coefficients should therefore 
be as small as possible, but still large enough to damp any instabilities. Although optimum values will vary 
from problem to problem, recommended levels are tf = 0(1) and e 7 = 2tf (Pulliam, 1986b). The recom- 
mended level for when used, is cf = 0(1). 

9.2 NONLINEAR COEFFICIENT ARTIFICIAL VISCOSITY 


The nonlinear coefficient artificial viscosity model is strictly explicit. Using the model as described by 
Pulliam (1986b), but in the current notation, the following terms are added to the right hand side of 
equation (8.5a). 


V 


f 



( 7 ) (+ , + ( 7 ),]('?’ M - 4VsV»,} 
( T ). +l + ( T - <VAQ)J 


(9.4) 


The difference operation A,V { A { Q is given by 


A i v f A jQ< “ Q/+ 2 ™ 3 Q<+i + 3 Q/ ~ Q/-i 


In the expression (9.4), \p is defined as 


t = 'l'x+ ty 


(9.5) 


10 These equations represent the use of the constant coefficient artificial viscosity model presented in this section. The 
nonlinear coefficient model to be presented in Section 9.2 is more complicated, but the same principle applies. 
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where \j/ x and \p y are spectral radii defined by" 

\U\+ay/tl+Z 2 y 

* x ~ (9.6) 

I v\ + aj t}l + n 2 y 

Here U and V are the contravariant velocities without metric normalization, defmed by 

U = + £x u "b %y v ^9 7) 

V = i/ r + r} x li + rjyV 


and a = JyRT , the speed of sound. 

The oarameters e® and e<«> arc the second- and fourth-order artificial viscosity coefficients. Instead of 
being specified directly by the user, as they are in the constant coefficient model, in the nonlinear coefficient 
model they « function of the pressure field. For the coefficients of the « direction differences, 

(9.8a) 


(e| 2) ). = k 2 At max(<r i+1 , <t,_i) 


= ma\[0, k 4 At - 


(9.8b) 


where 


Pi + 1 - 2 Pi+Pi-i 
p i+ 1 + 2 Pi+Pi~\ 


(9.9) 


Similar formulas are used for the coefficients of the jj direction differences. 

np« : s a nressure gradient scaling parameter that increases the amount of second-order 

JSE53S lUSS ™oo.hing neaj Lck wave. Tl,e logic used compute swuches 
off the fourth-order smoothing when the second-order smoothing term is large. 

The oarameters k, and k 4 are user-specified constants. Like the coefficients in the constant coefficient 
model the optimum values will be problem-dependent, and are best chosen through expenence. Cases have 
teen run wUhTadSs of k 2 ranging from from 0.01 for flows without shocks to 0.1 for flows with shocks 
and ic 4 ranging from 0.0002 for flows computed with spatially constant second-order tira differencmg 
0.005 for flows computed with spatially varying first-order time differencmg. Pulliam (1986b) gives 
kj = 0.25 and k 4 = 0.01 as typical values for an Euler analysis. 

Like the constant coefficient artificial viscosity model, the nonlinear coefficient model requires special 
formulas near boundaries. To apply (9.4) at i = 2, tf is needed at i = 1. It is defmed as 

(e^), = k 2 At max(<r 2 . <*i) 

With the above definition, applying (9.4) at i = 2 and / = N x - 1 requires * at / = 1 and i = N v They are 
defined as 


*?=£= T wtilc y in PROTEUS A« - 1/(N, - 1) and A, = 1/M 1)- The definitions used here for and 
result in an artificial viscosity level equivalent to that described by Pulliam. 
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~P± + 4 Pi - $P2 + 2^1 


Pa + ^P3 + 5p 2 + 2p l 

~? N \ -3 + 4 /V, -2 ~ -Va'i -i + 2/fy 
-3 + 4 Av 2 -2 + 5 /V, -i + 2/> V| 

And, finally, applying (9.4) at /= 2 and i = (V, - 1 requires AVAOnz-l a „,i •_ m 
numerous ormulas that could be used. I'he ones currently in tfieVko TliUS code ai ' 

A ^A,Q, = -Q 5 + 5Q 4 - 9Q 3 + 7Q 2 - 2Q, 

A S V i\Q/f t -i = Q,V, -4 - 5Q v, _3 + 9 Q*. -2 - 7Qa- + 2Q iVi 


are 


1 . 1 here are 
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APPENDIX A - EXPANSION OF VISCOUS TERMS 


In Section 5.2, the viscous terms in the governing equations arc linearized. To do this, the elements of 

E k and FV, given in equations (2. 1 7d) and (2. 1 7e) must first be rewritten in terms of the dependent variables, 
and with derivatives in the Cartesian directions transformed to derivatives in the computational directions 

using the chain rule. The non-cross derivative terms, involving E K and F„ , are then linearized using Taylor 

A A 1 1 

series expansion. The cross derivative terms, involving E y 2 and F are simply lagged one time level. This 
Appendix presents the fully expanded viscous terms required in the linearization procedure. 


The viscous term is given by equation (2.17d), which is repeated here. 


A 



1 

J 


1 

Re r 


0 

r xx^x T r xy^y 
T xy^x T T yy^y 

Pxtx + Pyty 


(A.l) 


where 

*** = 2 ^ u x + 2 (“x + v y) 

Tyy ~ T T Vy) 

*xy = + v*) 

P x = xx T VT xy p r 4x 

Py ~ UJ xy + vr yy p r Qy 

<7* = - kT x 
q y = —kT y 

The chain rule is used to transform derivatives in the Cartesian directions into derivatives in the com- 
putational directions, resulting in 

T xx = T T T d~ tj y^rf) 

~tyy = ( 2 /i T T d" *lx^r}) 

T xy d T d 

P x = (2 /x + K){Zx uu i + + + *?>%) 

+ m( V w ? + -I- r, x w n ) + — (^7 { + y] x T v ) 

P y = (2/i + 2)(£yW { + ijyW^) + + rixvuj 

Jc 

+ + VyUUr, + + 1^’,) + ^ (fj- 7 * + 
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The above expressions for the ts and are next substituted into equation (A.l). The f derivative terms 

A A 

become elements of E K , and the rj derivative terms become elements of E,, . The resulting four elements 
of E V[ (excluding the 1 /JRe, coefficient) are 

(E^)i = 0 (A. 2a) 

(E Ki ) 2 = 2 ntlu f + H x (£xU f + £ y Vj) + fity({ y u s + Z x v s ) (A. 2b) 

+ Xyit^ + ( y v t ) + + t x v s ) (A. 2c) 

(E k ,) 4 = 2 + ^w { ) + + ^uv { ) + 2^(^vu { + ^w { ) 

+ ntx(ty vu ! + ^X w {) + tfyityUUf + (A. 2d) 

For linearization it is convenient to rewrite the last element as 

(E„,)4 = -y-- iiW)( + t 2 y( v \l + (#* + 

+ y [&v 2 ) { + «J(« 2 )j] + yr (£ + (A.2e) 

A A 

The elements of F Kj have exactly the same form as those of E but with £ replaced by rj . 

The four elements of Ey 2 (again excluding the 1 jJRe r coefficient) are 

(E^ 2 )i = 0 (A. 3a) 

A 

(^^2)2 x^r] *F (Wx^rj "F "F ^^yi^ly^rj "F Wx^rj) (A. 3b) 

A 

(£^2)3 T ^£y{Vx^t] *F ^ly^r)) "F ^^x^y^rj ~F *7x v >/) (A. 3c) 

A 

(Ej /^)4 2fj.{^ x rf x uu^ 4- ^ytjyWy) 4- X^ x {yj x uu^ 4~ 4- A^y{yj x vu^ 4~ 

k. 

+ nZJfly™* + 1x w *) + nt/'lyWr, + + y (ijc»Jx + (A. 3d) 

1 r r 

The elements of F K have exactly the same form as those of E„ , but with £ replaced by r\ and yj replaced 
by <T 
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APPENDIX B - AXISYMMETRIC ANALYSIS 


The analysis usd PROTljUS for rip** 

r,£ r P Sc y C? ssl, rw «* ««; -« .*>— -** is 

described separately in this appendix. 


B.l GOVERNING EQUATIONS 

In cylindrical coordinates, the governing equations for axisymmetric flow, with swirl, can be written 
using vector notation as 

0 (rQ) 3 (rE j j _ + 3 (rF^ + H r (B.l) 

dt + dx dr dx dr 


where 


Q = [p pu pv pw £ r ] 


E = 


pu 

pu +p 

puv 

puw 

(E T + p)u 


pv 

puv 

pv 2 +p 

pvw 

(E T +p)v 


H 


0 

0 

-p - p* 2 

pvw 

0 


(B. 2 a) 


(B.2b) 


(B.2c) 


(B. 2 d) 


1 XX 




‘xr 

T x0 

UT xx + VT xr + Wr x9 


Pr r 


<1x 


(B. 2 c) 
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0 


(B.2f) 


^xr + ^ rr + wr r0 -~ qr 


The shear stresses and heat fluxes are given by 


/-!« +±i 
dx r ' 

(^)] 

r -^-+±1 

dx r y 

d(rv) \1 

* )\ 

*L + ±{ 

d{rv) \"j 

dx + r ^ 

~)\ 

/ du dv 

\ 

V dr dx 4 

) 

dw 



r xO = H 


_ ( dyv w \ 


*--*■ f 

„ *” l ^ ese equatlons ’ JC ’ r ' , and0 represent the axial, radial, and circumferential directions respectively and 

jsssr The — *• *— « «* « 

For turbulent flow, M , 2, and * represent effective coefficients. The turbulence model is described in 

?;?• Jhe only modification to the model for axisymmetric flow is the definition of Ini the mag- 
nitude of the total vorticity. For axisymmetric flow, ' * ’ ag 
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When the generalized grid transformation of Section 2.3 (with y replaced by r), is applied to equation 
(B.l) the result is 

(r Q) t + (r Q)^ f + (r Q)^;, + (r ¥)rZ x + (r + (r F)^ r + (r F)^»j r + H 

- (r E v ) f i x - ( r E v ) n n x - (r F„) { {, - (r Fy)^, — Hy— 0 (B-4) 

Although this axisymmctric equation cannot be put into exact strong conservation law form, the pro- 
cedure used to do so for the two-dimensional equation, described m Section 2.4, is nonetheless applied to 
equation (B.4). ihe result is 


d(r Q) d(r E) d(r F) a _ djr E^) d(r¥y) £ 

~a7~ + a* + a,, + at dr, 1 


(B.5) 


where 


Q = 


Q 

J 


E = -t-(E* x + F<;, + Q£,) 

F = -y (Ei/ X + Fi? r + Qrj,) 

H = f 

E[/= -y (V.yZ x + ¥ V %r) 

¥y = -j(V.y*, x + Fyt, r ) 

A H y 

Ho- 
using equations (B.2a) through (B.2g) these can be expanded as 

Q = -y [p pu pv pw E T ] r 


E 


1_ 

J 


pu£ x + - pvZ r + PZ, 

(pu 2 + p)Z x + + P u $t 

pM'i x + (P v2 + p)i r + pv£ t 
puw^x + pvw£ r + pw£, 

( E r + p) u£ x + (E r + p)v£ r + E t £, 


A 

F 


pur\ x + pvt , r + pr\, 

(pu + p)t] x + puvt , r + pur,, 
puvr, x + (pv 2 +p)r, r + pvt,, 
puwt\ x + pvwt, r + pwtjf 
(E T + p)ut, x + (E T + p)vr, r + E T t\ t 


(B.6a) 


(B.6b) 


(B.6c) 
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11 = 


0 

0 

-p - pw 2 

pvw 

0 


(B.6d) 


A 

K. 


1 1 
J Re r 


1 1 
J Re r 



0 


T xx$x 

+ 

T xrZr 

T xrZx 

4- 

?rrZr 

r x0% x 

4- 

T r0^r 

Pxix 

+ 

Prir 


0 


X xx^x 

4- 

T xr y lr 

r xr>,x 

+ 

T r r r,r 

T x0*1x 

4- 

r rO y lr 

PxVx 

4- 

PrVr 


(B.6c) 


(B.6f) 


where 


A 1 I 


0 

0 

~ T oo 

T rO 

0 


(B.6g) 


fix = UT XX + -y^Vx 

Pr = UT xr + "I rr + ~ J^T <7r 


B.2 LINEARIZATION 


(B.7) 


A 

Solving equation (B.5) for dQ/dr (assuming r is not a function of time) and substituting the result into the 
time differencing scheme of Beam and Warming, given by equation (4.1), for d{A(y')jd-c and dQ'jd r yields 

_+i+i if 5^L> I I A ,yA Aj_xf | ,y\ 

i + o 2 ' ( d( a, i J i + o 2 ' ( d( + a, + 11 ) 


0 , At if d(rAK n y ) d(rAl-'l) 




+ 


1 +0 2 r 

0 


di 


dr, 


1 + 0 , r 


dZ, 


dr, 


1+0 


2 AC/ 2-1 + O 


2 


( 0 ,-y- 0 2 )(Ar) 2 + (AT) 3 J 


(B.8) 


This equation must be linearized using the procedure described in Section 5.0. 
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B.2.1 Inviscid Terms 


Tor the inviscid terms the Jacobian coefficient matrix <3E/dQ is 


dQ 


~ w f\ 


?,+/ i+“{* + 


dp 

d(pu) 




i, 

dp . 
+ ao»v) 






V js 4- — l 

** a(p«) 


-4-t) 


/ 2 {.+/l 


h d(pv) 

d(pw) 


t,+A 

dp 

f dp 

d(pv) 

d(pw) 


a£ r 

ap 

d£ r 




(B.9) 


where y; = u{ x + v£, and f 2 = (E T + p)lp. The Jacobian matrix dF/dQ has the same form as dE/dQ, but with 
£ replaced by »j. 

For the additional term H, the linearization procedure gives 


an 

A 

dQ 


0 

0 

^ 2 

dp 

—vw 


0 

dp 

d(pu) 

0 

0 


0 

dp 

d(pv) 

w 

0 


0 


dp 


d(pw) 

v 


— 2w — 


0 


dp 

dE r 

0 

0 


(BIO) 


B.2.2 Viscous Terms 

To linearize the viscous terms, F. y , K y , etc., must first be rewritten in terms of the dependent variables, 
and with derivatives in the cylindrical coordinate directions transformed to derivatives in the computational 
directions using the chain rule. The shear stress and heat flux terms, given by equations (B.3) and (B.7), 
become 

Tjw = (2 n + + Vx“r,) + T CW")e + 

T rr = 2u(b r v > + r lr V rj) + + *Jjc U rj) + ~f 

r e0 = 2p Y + >!(£,«{ + r\ x u n ) + T [W"){ + hi”)*] 

T X r = p(Z r u S + *lr u r, + £jc v { + Ix^) 

TxB = + ^x w V ) 
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T rl) = ti(Zr w S + >lr w n ) ~ b ~J~ 


Px = ( 2 1* + 2)(C x uu Z + 1x uu r,) + 4 [^ w ( n ’)j + >7r^(^), 7 ] 

+ + t? r v^ + £ X W, + + ^mi^) 

+ -j^(txT t + v x TJ 

Pr = + »/r vv n ) + "4 U>'M 7 + *1r v ( n ')r,] 

+ H(ir ltu > + V% + ^x w l + } 1x ltv n ) + l l (*r ww i + >lr ww v ) 

+ + >7x%) -b^ + -~ (trT s + nr'Q 

The above expressions for the shear stress and heat flux terms are substituted into equations (B.6e) 
through (B.6g). As in the two-dimensional planar case, the cross derivative terms are separated from the 
non-cross derivative terms. In addition, for the axisymmetric case the non-derivative terms are included 
with the cross derivatives. 


The resulting five elements of E v (excluding the \jJRe r coefficient) are 

(tV,)t = o 

(kfqb = 2/r£xW* + + T *r( n ')|] + t l *r(£r u r + £jc v j) 

( E ^,)3 = + 4- + bixdr u i + Cjd’f) 

(E [/,).! = m£x w { + b^r w r\ 

( E ^,)5 = 2m(£* mu { + + H x £x uu > + “4 + 2^ r ^x VM f + ~Y £r v ( n ’){l 

+ H$ x (t r VU { + + ^MW { ) + + ( x uv g + t r WW g ) + -£-(£+ 


Pr r 

For linearization it is convenient to rewrite the last clement as 

= — 2 — >e + & v ){] + (M + *K x tr(wk + K ~r (& + **«v) 

+ y + ” 2 ) { + Z?(“ 2 + * 2 ) f ] + -^T (d + d)7> 

I he elements of FV, have exactly the same form as those of E K , but with £ replaced by >/. 
The five elements of E y 2 (again excluding the 1 \JRe, coefficient) are 

(E^), =0 

i^V 2 h = ^x*lx u n + 2^x[ > fx w r, + ~r VrM,,] + ^ri.Vr u n + Vx v n ) 


(B.lla) 

(B.llb) 

(B.llc) 

(B.lld) 

(B.lle) 


(B.llf) 


(B. 12a) 
(B. 12b) 
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(Ek 2 )3 = 2 l l Zr>1r v n + Hr[yx u r, + T 'b^J + l l *x(lr“ n + >lx v r, ,) (B. 1 2c) 

(E(/ 2 )4 = l l ix r lx w rj + tfr1r w n ~ "7" (B.12d) 

(E^S = M^x*lx uu n + trtr w J + + T J + + T frK") J 

+ ^x(Vr vu v + 1x w r, + »?**%) + vtri'h 1 ^ + ’/jc^ + >//■*%) 

- -V" + (£jc*1jc + ZrlrWr, (B.12c) 

The last element can he rewritten as 

(E|/ 2 )5 = MtxVxWr, + ttfr^r,) + 2 Zx('1x UU r l + 1r w „) + / ‘^r( f lx vu n + >»,%) + *'b T ^x u + Zr v ) r v 

+ n* x (>1r v u r, + *»x W n + >?x M ' vv f ,) + H^ r {n r uu n + »7*% + '»#•*%) 

2 

-Mtr^r + -£r(tx , lx+ $r'1r) T r t (B.121) 

A A 

The elements of F v have exactly the same form as those of E Kj , but with C replaced by >/ and >/ replaced 
by £. 

The five elements of H v are 

(H„), =0 (B. 13a) 

(lV )2 = ° (HA3b) 

(H„) 3 = -2fi y - + ri x u n ) + — + tirirv)^ (B. 13c) 

(Hi) 4 = M€ r w { + *i r w J -fy (B. 13d) 

(lV)s = 0 (B.13e) 
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Performing the linearization, the Jacobian coefficient matrix dY^JdQ is 


where 


dF. v 


dQ 


1 

Re r 



0 



0 

0 

0 


3 ( 

» 

^ t 

a "ir\ 

7 ) f *'« 7- f < 

0 

0 


3 ( 

» 

•-i< 

,b) + “" T r i 

0 

0 

(B 14 ) 



*xx = (2/i + 

a rr - PZx + ( 2 / 4 + 


= "b M^r 
a xr = (M “b 
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<*v, 


at: 


E, \ , 5 

+ «0 


(3Q y 52 \ dQ }i\ 

BE,, \ ( d El 




(-M-) 

V */>«) / 


<?Q /53 


3E 


Ei 


dQ / 31 


3E t 


+ a rr "p" A 


+ «0 


. d ( 3T \ 
{ + «o ^ ^ d{pv) ) 


dQ /54 


<3Q / 4 1 




V 3(pw) ) 


The Jacobian coefficient matrix for the remaining non-cross derivative viscous terms, dlyJdQ, has the 
same form as dE v JdQ, but with Z replaced by tj. 

And finally, linearizing H v , the Jacobian coefficient matrix dH K /dQ is 


where 


dn v 

A 

dQ 


Re, 


5H 


v 


dQ 

dHy 


31 


dQ U\ 


0 

dHy 

dQ / 32 

0 


dHi 

^5 733 


aH, 


o 


dQ J 44 
0 0 


an ( 

dQ J 3i 


*ix 


d 

dt 




+ [2p + m r r f 'f ’iSr }) 3 r p ^ r 3»j ( P ) 


d H t 

<3Q / 32 




_a_ 
=JC a£ 


(^)33 = _ ^ r ^(^)~ C2M + ;( ^ + nA)] ^^ _ ^^(^^ 


(R 15) 
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OQ ) 41 


E d / w \ I 1 w d / w \ 

^'dF\T) + -Fir-^r^-( -p) 


(li v \ _ * JL/ J_\ a 1 . a / 

<30 / 44 * V P ) rP+^'dA 


B.2.3 liquation Of Stale 


The equation of state given in Section 5.3 must be modified slightly to add the swirl velocity vv. Thus, 


P — (y ~ 1) 2 p{ u ^ t v " + ;y2 )J 


or, in terms of temperature, 


r.-lTiz: 
<v L p 


1 / 2 . 2 , 2 \ 
" 2 " (m + v -t- vv J 


The derivatives arising from the linearization are the same as those presented in Section 5.3, except for 


dp _ y - \ 
dp 2 


(« 2 + v 2 + w 2 ) 


= -(>■- l)w 


dT _ 1_ JV 1 / 2 2 2- 

T — r 1 7T ( ^ Tv T VV 

Op Cy l p \ 


(B.18c) 


d(pw) C V P 

11 constant stagnation enthalpy can be assumed, the appropriate equation of state is 


P--y-p 


[ A /-y( w 


2 + V 2 + vv 2 ) J 


and the temperature becomes 


1 '= [/‘y - y (w 2 + V 2 + w 2 )J 


Again, the derivatives arising from the linearization are the same as in Section 5.3, except for 

d P v - 1 f, . 1 , 2 . 2 . 2.1 


= L -y— [a 7 + y (u 2 + v 2 + w 2 ) J 


dp y — 1 


(B.18d) 


(B.19) 


(B.20) 


( B . 2 1 a) 


(B.2 lb) 


d! 1,2.2, 2, 

— +v +w) 

dj_ = w 

d(pw) c pP 


(B.2 lc) 


(B.2 Id) 
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B.2.4 Linearized Governing Equation 

The linearized form of equation (B.8) can now be written as 


a 0, At } 

AQ" + 


1 +0 2 r 


d¥. 

A 

dQ 


AQ 


+ 


d n 


dF I A 

r l —r- 1 AQ 


5Q 


+ 


0|At 1 

s . 

ri 

/ A 

1 + 0 2 r < 

V ( 

\ 

v dQ 

At 1 ( 

✓ A 
f d(r E 

U 

d(rh 


1 +0 7 r 






AQ" 


+ H + 


dr] 


d¥ v Y .. 

y l I A 


A 

dQ 


AQ 


A 

<3H 

A 

<?Q 


+ 


AQ" 


dH, 


dQ 


AQ" 


At 1 ( d ( rE V? 5 ( rF ^,) 


dt 


dr, 


+ Hr 


(1 + 0 3 )At i / ^E t ' 2 ) Krty) \ 0 3 At i ( 3(rF^) 


t n — 1 


1 T $9 




+ 


dr] 


d( 


drj 


0 9 


AQ" -1 + O 


(fl, - y - ^ 2 )(At) 2 , (0 3 - 0i)(At) 2 , (At) 3 


l ~h O2 

BJ SOLUTION PROCEDURE 

Letting LHS(B.22) represent the left hand side of equation (B.22), we can write 


(B.22) 


LHS(B. 22 ) = < I + 


L 

1 + 0 2 ' 


AU_ r ^ 4 4-^ M an 

a { \ an 30 I dr] \ dQ dQ I 


d\\ v 


AQ (B.23) 


dQ dQ J \ dQ dQ I \ dQ dQ 

where I represents the identity matrix. The term in braces in equation (B.23) can be factored to give 


LHS(B. 22 ) = 


M J 1 d [ r gE 


a E^ 


n 


0 \At ) ( du d\\ v 


1 + 


1 + 0 2 r dl y d Q d Q 

d i A*z 1 d f dF _ d * y \ 


i + 0 2 r 


dQ dQ 


] + 0 2 r \ dQ dQ 


AQ 


( O t ST 

y . 

d 

( r jL 

A 

5 E„, 

\_a_| 

' r jL 

A 

aF ri 

VI 

V 1 + 0 2 . 

)v 

I ^ 

|<\> 

i 

( dQ 

“ r A 

dQ 

) dr, 1 

( aQ 

r A 

dQ 

)\ 


AQ 


( 

V 1 

( ILL 

dHy ^ 

|ii 

( * 

dF 

V 

V 1 + o 2 

)v 

V dQ 

I 

o>| 

) dr, | 

1 r A 

^ dQ 

dQ J 
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The last tw r o terms represent the splitting error. 

Equation (B.22) can thus be rewritten in spatially factored form, and, neglecting the temporal truncation 
and splitting error terms, becomes 
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Equation (B.25) can be split into the following two-sweep sequence 
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Applying the spatial differencing formulas of Section 6.0 results in 
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Sweep 2 (n direction) 
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These equations are solved using the same matrix inversion procedure described in Section 8.2. 
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